Loading required package: carData
Attaching package: 'car'
The following object is masked from 'package:psych':
logit
The following object is masked from 'package:dplyr':
recode
The following object is masked from 'package:purrr':
some
Code
# Read .csv file into a tibbledata <-read_csv("data/middle_school_music_mindsets_clean_partial.csv")
Rows: 503 Columns: 46
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (8): school, race, hispanic_latino, gender, currentSchoolMusic, school_...
dbl (38): survey_id, grade, schoolChoir, schoolInstrumental, schoolExtracurr...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Code
#reclassify character variables as factorsreclassFactor <-c("grade","school","race","hispanic_latino","gender","currentSchoolMusic","familyMusic")data[reclassFactor] <-lapply(data[reclassFactor], as.factor)
Code
# Reorder the levels for currentSchoolMusic and set "No music" as the reference leveldata$currentSchoolMusic <-relevel(data$currentSchoolMusic, ref ="No music")
I explored the predictors’ variance using Variance Inflation Factors (VIF) to identify potential collinearity. I followed the same steps for the three models (three different outcome variables).
Global MSC as outcome - VIF (Variance Inflation Factors)
Based on these results, musical_training, currentSchoolMusic, and schoolInstrumental should not be included in the model together due to their potential multicollinearity. To solve this issue, I removed one by one starting with the one with highest VIF, currentSchoolMusic.
After removing currentSchoolMusic and musical_training, there seems to be no collinearity between the predictors. Hence,I will not include those variables in the model.
Ability MSC as coutcome - VIF (Variance Inflation Factors)
Similar to the previous model, musical_training, currentSchoolMusic, and schoolInstrumental should not be included in the model together due to their potential multicollinearity. To solve this issue, I removed one by one starting with the one with highest VIF, currentSchoolMusic.
After removing currentSchoolMusic and musical_training, there seems to be no collinearity between the predictors. Hence,I will not include those variables in the model.
Non-Academic MSC as Outcome - VIF (Variance Inflation Factors)
Consistent with the first two models, musical_training, currentSchoolMusic, and schoolInstrumental should not be included in the model together due to their potential multicollinearity. To solve this issue, I removed one by one starting with the one with highest VIF, currentSchoolMusic.
After removing currentSchoolMusic and musical_training, there seem to be no collinearity between the predictors. Hence,I will not include those variables in the model.
Assumptions
Before building the linear regression models, I examined assumptions of linearity, homoscedasticity, and other assumptions. I also examine if using the music self-concept composite score as an outcome would yield different results compared to the musical ability/music self-concept subscale.
Kolmogorov-Smirnov and Shapiro- Wilk tests showed that the data follows a nonnormal distribution. Individual examination of each variables’ skewnewss and kurtosis identify three variables, but the distribution was not extreme.
Stepwise Multiple Linear Regression
I conducted a series of linear regression models to examine the relationships between multiple predictors and each of the three outcome variables, global MSC, ability MSC, and non-academic MSC. The aim was to address the research questions gradually, facilitating the interpretability of the findings by assessing the influence of successive categories of predictors on the models’ improvement. I conducted the same steps and used the same predictors for each outcome variable.
On Step 1, I introduced mindset of music ability as the sole predictor.
To investigate the impact of various types of music engagement on music self-concept, I introduced two sets of variables in the following steps. Step 2 included active engagement, whereas Step 3 comprised variables related to duration of engagement in music activities in and out of school.
In the following and final step (Step 4), I introduced background variables like grade, school, race, ethnicity, gender, and family music background.
GLOBAL MSC - Stepwise Multiple Linear Regression
Intercept
Code
#Interceptlm_msc_full_step_0 <-lm(formula = msc_full ~1, data = data)summary(lm_msc_full_step_0)
Call:
lm(formula = msc_full ~ 1, data = data)
Residuals:
Min 1Q Median 3Q Max
-44.748 -8.748 0.252 9.252 32.252
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 75.7475 0.6215 121.9 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.94 on 502 degrees of freedom
Call:
lm(formula = msc_full ~ growth_mindset + active_engagement, data = data)
Residuals:
Min 1Q Median 3Q Max
-36.328 -5.335 0.328 5.410 26.347
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.31352 2.20194 7.863 2.33e-14 ***
growth_mindset 1.79947 0.14887 12.087 < 2e-16 ***
active_engagement 1.56372 0.07537 20.748 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.752 on 500 degrees of freedom
Multiple R-squared: 0.6073, Adjusted R-squared: 0.6057
F-statistic: 386.6 on 2 and 500 DF, p-value: < 2.2e-16
Code
# Calculate model fit improvementf_test_result_lm_msc_full_1_2 <-anova(lm_msc_full_step_1, lm_msc_full_step_2)f_test_result_lm_msc_full_1_2
Analysis of Variance Table
Model 1: msc_full ~ growth_mindset
Model 2: msc_full ~ growth_mindset + active_engagement
Res.Df RSS Df Sum of Sq F Pr(>F)
1 501 71272
2 500 38299 1 32974 430.48 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
In step 2, active engagement emerged as a significant predictor, collectively contributing to 62% of the variance in MSC (adjusted R2 = 0.6057, p < .001 ). The addition of these variables significantly improved the model, F = 430.48, p < .001).
The refined model explained 65% of the variance in music self-concept (adjusted R2 = 0.6404, p < .01), representing a significant fit improvement, F = 7.89, p < .001).
In the new model, gender was a significant predictor of global music self-concept, with male students differing significantly from female students (reference group). This extended model explained 68% of the variation in music self-concept (adjusted R2 = 0.6455, p = .049), showing a significantly superior fit compared to the previous model (F = 1.55, p < .01).
Results Table
Code
#RESULTS TABLE Global MSC# Extract coefficients and statisticsresults_lm_msc_full <-bind_rows(broom::tidy(lm.beta(lm_msc_full_step_4), conf.int =TRUE))# Round output to two decimal pointsresults_lm_msc_full$estimate <-round(results_lm_msc_full$estimate,2)results_lm_msc_full$std_estimate <-round(results_lm_msc_full$std_estimate,3)results_lm_msc_full$std.error <-round(results_lm_msc_full$std.error,2)results_lm_msc_full$statistic <-round(results_lm_msc_full$statistic,2)results_lm_msc_full$p.value <-round(results_lm_msc_full$p.value,3)results_lm_msc_full$conf.low <-round(results_lm_msc_full$conf.low,2)results_lm_msc_full$conf.high <-round(results_lm_msc_full$conf.high,2)# Combine conf.low and conf.high into a single column called "95% CI" and change column namesresults_lm_msc_full <- results_lm_msc_full %>%mutate("95% CI"=paste0("[", conf.low, ", ", conf.high, "]")) %>%select(-conf.low, -conf.high) %>%# Remove the original conf.low and conf.high columnsrename(b = estimate, β = std_estimate,SE = std.error,t = statistic,p = p.value)# Change order of columnsresults_lm_msc_full <- results_lm_msc_full %>%select(term, b, SE, β, t, "95% CI", p)# Change output for p values undeer .001results_lm_msc_full <- results_lm_msc_full %>%mutate(p =ifelse(p <0.001, "<.001", p))# Change predictors' namesresults_lm_msc_full <- results_lm_msc_full %>%mutate(term =case_match(term, "(Intercept)"~"(Intercept)","growth_mindset"~"Growth mindset", "active_engagement"~"Active Engagement","schoolChoir"~"School choir","schoolInstrumental"~"School instrumental","schoolExtracurricular"~"School music club","outsideChoir"~"Choir (outside school)","outsideInstrumental"~"instrumental (outside school)","privateLesson"~"Private lesson","selfTaught"~"Self-teaching","grade7"~"Grade 7","grade8"~"Grade 8","schoolB"~"School B","schoolC"~"School C","familyMusicYes"~"Family music engagement","raceBlack or African American"~"Black or African American","raceSome Other Race"~"Other race","raceTwo or More Races"~"Two or more races","raceWhite"~"White","hispanic_latinoYes"~"Ethnicity","genderMale"~"Male","genderNon-Binary"~"Non-Binary","genderPrefer not to say"~"Gender not disclosed") )print(results_lm_msc_full, n=100)
# A tibble: 23 × 7
term b SE β t `95% CI` p
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 (Intercept) 24.4 3.2 NA 7.62 [NA, NA] <.001
2 Growth mindset 1.58 0.16 0.309 9.87 [-0.01, 0.62] <.001
3 Active Engagement 1.3 0.08 0.504 15.7 [0.34, 0.67] <.001
4 School choir 0.15 0.23 0.02 0.66 [-0.43, 0.47] 0.512
5 School instrumental -0.2 0.26 -0.023 -0.78 [-0.53, 0.48] 0.433
6 School music club 0.58 0.29 0.062 2 [-0.51, 0.63] 0.047
7 Choir (outside school) 0.44 0.43 0.03 1.03 [-0.82, 0.88] 0.305
8 instrumental (outside school) -0.31 0.49 -0.019 -0.64 [-0.99, 0.95] 0.524
9 Private lesson -0.07 0.21 -0.009 -0.31 [-0.43, 0.41] 0.757
10 Self-teaching 1.26 0.23 0.176 5.43 [-0.28, 0.63] <.001
11 Grade 7 -1.16 0.97 -0.039 -1.2 [-1.94, 1.86] 0.231
12 Grade 8 -1.42 0.97 -0.049 -1.47 [-1.95, 1.85] 0.142
13 School B -0.72 0.95 -0.026 -0.76 [-1.89, 1.84] 0.45
14 School C 0.4 1.07 0.013 0.37 [-2.1, 2.12] 0.711
15 Family music engagement -0.08 0.78 -0.003 -0.11 [-1.54, 1.53] 0.916
16 Black or African American 2.52 2.43 0.05 1.04 [-4.72, 4.82] 0.299
17 Other race 1.09 2.77 0.015 0.39 [-5.43, 5.46] 0.694
18 Two or more races 1.24 2.52 0.023 0.49 [-4.92, 4.97] 0.622
19 White -0.21 2.06 -0.006 -0.1 [-4.06, 4.05] 0.919
20 Ethnicity -0.52 1.89 -0.008 -0.27 [-3.72, 3.7] 0.784
21 Male -2.48 0.82 -0.089 -3.03 [-1.69, 1.52] 0.003
22 Non-Binary -2.95 4.28 -0.019 -0.69 [-8.43, 8.39] 0.491
23 Gender not disclosed -5.9 3.03 -0.053 -1.95 [-6.01, 5.9] 0.052
The overall significance of the final model was confirmed through the F-test (F = 42.55, p < .01).
Post Hoc Analyses
Equality of coefficients (Wald Test)
Although several predictors showed significant effect on the outcome, their coefficient may or may not differ from each other significantly. To examine if the differences between the coefficients (as observed in values of b and β ) are significant, I visually inspected those coefficients.
Visual inspection of the data suggests coefficients for mindset of music ability, active engagement, and self-teaching may not be significantly different. Hence, I investigated further by conducting Wald Tests between each pair of predictors.
Linear hypothesis test:
active_engagement - selfTaught = 0
Model 1: restricted model
Model 2: msc_full ~ growth_mindset + active_engagement + selfTaught +
gender
Res.Df RSS Df Sum of Sq F Pr(>F)
1 497 34290
2 496 34276 1 13.72 0.1985 0.6561
Wald tests results indicated that the difference in coefficients between mindset of music ability and active engagement, mindset of music abilty and self-teaching, and active engagement and self-teaching were not statistically significant, suggesting those predictors had the same effect on students’ global music self-concept.
Joint Significance (Wald Test)
Code
wt_gender <- lmtest::waldtest(lm_full, term ="gender")print(wt_gender)
I employed the Wald test once more to delve deeper into the joint significance of the gender variable. The results showed a significant difference between the model incorporating gender and the one excluding it, F = 4.2502945, p = < .01. This finding indicates that gender stands as a significant predictor of music self-concept, underscoring its importance in the predictive model.
Scatter Plots
Code
#scatter plots with fitted regression linesplot1_full <-ggplot(data, aes(x = growth_mindset_z, y = msc_full_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Mindset of Music Ability",y ="Overall") +theme_ipsum()plot2_full <-ggplot(data, aes(x = active_engagement_z, y = msc_full_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Active Engagement",y ="Global Music Self-Concept") +theme_ipsum()plot3_full <-ggplot(data, aes(x = selfTaught_z, y = msc_full_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Self-Teaching",y ="Global Music Self-Concept") +theme_ipsum()plot4_full <-ggplot(data, aes(x = schoolExtracurricular_z, y = msc_full_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Extracurricular School Music Clubs",y ="Global Music Self-Concept") +theme_ipsum()#Combine plots into a single display (grid layout)grid.arrange(plot1_full, plot2_full, plot3_full, plot4_full, ncol =2)
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
ACADEMIC MSC (Music Ability Only) - Stepwise Multiple Linear Regression
Intercept
Code
#Interceptlm_msc_ability_step_0 <-lm(formula = musical_ability_msc ~1, data = data)summary(lm_msc_ability_step_0)
Call:
lm(formula = musical_ability_msc ~ 1, data = data)
Residuals:
Min 1Q Median 3Q Max
-7.173 -3.173 -0.173 2.827 7.827
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.1730 0.1603 75.92 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.596 on 502 degrees of freedom
Call:
lm(formula = musical_ability_msc ~ growth_mindset + active_engagement,
data = data)
Residuals:
Min 1Q Median 3Q Max
-7.5612 -1.8720 0.0405 1.8709 9.4115
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.37141 0.69821 -0.532 0.595
growth_mindset 0.66545 0.04721 14.097 < 2e-16 ***
active_engagement 0.18323 0.02390 7.667 9.24e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.775 on 500 degrees of freedom
Multiple R-squared: 0.4068, Adjusted R-squared: 0.4045
F-statistic: 171.5 on 2 and 500 DF, p-value: < 2.2e-16
Code
# Calculate model fit improvementf_test_result_lm_msc_ability_1_2 <-anova(lm_msc_ability_step_1, lm_msc_ability_step_2)f_test_result_lm_msc_ability_1_2
Analysis of Variance Table
Model 1: musical_ability_msc ~ growth_mindset
Model 2: musical_ability_msc ~ growth_mindset + active_engagement
Res.Df RSS Df Sum of Sq F Pr(>F)
1 501 4303.4
2 500 3850.7 1 452.72 58.784 9.239e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculate model fit improvementf_test_result_lm_msc_ability_2_3 <-anova(lm_msc_ability_step_2, lm_msc_ability_step_3)f_test_result_lm_msc_ability_2_3
# Calculate model fit improvementf_test_result_lm_msc_ability_3_4 <-anova(lm_msc_ability_step_3, lm_msc_ability_step_4)f_test_result_lm_msc_ability_3_4
#RESULTS TABLE: Ability MSC# Extract coefficients and statisticsresults_lm_msc_ability <-bind_rows(broom::tidy(lm.beta(lm_msc_ability_step_4), conf.int =TRUE))# Round output to two decimal pointsresults_lm_msc_ability$estimate <-round(results_lm_msc_ability$estimate,2)results_lm_msc_ability$std_estimate <-round(results_lm_msc_ability$std_estimate,2)results_lm_msc_ability$std.error <-round(results_lm_msc_ability$std.error,2)results_lm_msc_ability$statistic <-round(results_lm_msc_ability$statistic,2)results_lm_msc_ability$p.value <-round(results_lm_msc_ability$p.value,3)results_lm_msc_ability$conf.low <-round(results_lm_msc_ability$conf.low,2)results_lm_msc_ability$conf.high <-round(results_lm_msc_ability$conf.high,2)# Combine conf.low and conf.high into a single column called "95% CI" and change column namesresults_lm_msc_ability <- results_lm_msc_ability %>%mutate("95% CI"=paste0("[", conf.low, ", ", conf.high, "]")) %>%select(-conf.low, -conf.high) %>%# Remove the original conf.low and conf.high columnsrename(b = estimate, β = std_estimate,SE = std.error,t = statistic,p = p.value)# Change order of columnsresults_lm_msc_ability <- results_lm_msc_ability %>%select(term, b, SE, β, t, "95% CI", p)# Change output for p values undeer .001results_lm_msc_ability <- results_lm_msc_ability %>%mutate(p =ifelse(p <0.001, "<.001", p))# Change predictors' namesresults_lm_msc_ability <- results_lm_msc_ability %>%mutate(term =case_match(term, "(Intercept)"~"(Intercept)","growth_mindset"~"Growth mindset", "active_engagement"~"Active Engagement","schoolChoir"~"School choir","schoolInstrumental"~"School instrumental","schoolExtracurricular"~"School music club","outsideChoir"~"Choir (outside school)","outsideInstrumental"~"instrumental (outside school)","privateLesson"~"Private lesson","selfTaught"~"Self-teaching","grade7"~"Grade 7","grade8"~"Grade 8","schoolB"~"School B","schoolC"~"School C","familyMusicYes"~"Family music engagement","raceBlack or African American"~"Black or African American","raceSome Other Race"~"Other race","raceTwo or More Races"~"Two or more races","raceWhite"~"White","hispanic_latinoYes"~"Ethnicity","genderMale"~"Male","genderNon-Binary"~"Non-Binary","genderPrefer not to say"~"Gender not disclosed") )print(results_lm_msc_ability, n=100)
# A tibble: 23 × 7
term b SE β t `95% CI` p
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 (Intercept) 1.98 0.91 NA 2.18 [NA, NA] 0.03
2 Growth mindset 0.41 0.05 0.31 8.98 [0.22, 0.4] <.001
3 Active Engagement 0.13 0.02 0.19 5.41 [0.15, 0.24] <.001
4 School choir 0.09 0.07 0.04 1.36 [-0.08, 0.17] 0.175
5 School instrumental 0.3 0.07 0.13 4.04 [-0.01, 0.28] <.001
6 School music club 0.37 0.08 0.15 4.53 [-0.01, 0.32] <.001
7 Choir (outside school) 0.19 0.12 0.05 1.53 [-0.19, 0.29] 0.126
8 instrumental (outside school) -0.41 0.14 -0.1 -2.97 [-0.37, 0.18] 0.003
9 Private lesson 0.19 0.06 0.11 3.2 [-0.01, 0.23] 0.001
10 Self-teaching 0.44 0.07 0.24 6.72 [0.11, 0.37] <.001
11 Grade 7 -0.55 0.27 -0.07 -1.99 [-0.61, 0.47] 0.047
12 Grade 8 -0.95 0.27 -0.13 -3.45 [-0.67, 0.41] 0.001
13 School B -0.28 0.27 -0.04 -1.03 [-0.57, 0.49] 0.303
14 School C 0.85 0.3 0.1 2.79 [-0.49, 0.7] 0.005
15 Family music engagement 0.28 0.22 0.04 1.25 [-0.4, 0.47] 0.211
16 Black or African American 0.18 0.69 0.01 0.26 [-1.34, 1.37] 0.796
17 Other race -0.25 0.79 -0.01 -0.32 [-1.56, 1.53] 0.75
18 Two or more races -0.02 0.71 0 -0.03 [-1.4, 1.4] 0.974
19 White -0.68 0.59 -0.08 -1.16 [-1.23, 1.07] 0.247
20 Ethnicity -0.91 0.54 -0.05 -1.69 [-1.11, 1] 0.091
21 Male 0.71 0.23 0.1 3.07 [-0.36, 0.55] 0.002
22 Non-Binary 0.55 1.21 0.01 0.45 [-2.37, 2.4] 0.65
23 Gender not disclosed -1.45 0.86 -0.05 -1.69 [-1.74, 1.64] 0.092
# Calculate model fit improvementf_test_result_lm_msc_ability_4_interactions <-anova(lm_msc_ability_step_4, lm_msc_ability_interactions)f_test_result_lm_msc_ability_4_interactions
# Calculate model fit improvementf_test_result_lm_msc_ability_4_interactions_short <-anova(lm_msc_ability_step_4, lm_msc_ability_interactions_short)f_test_result_lm_msc_ability_4_interactions_short
#Plotinteract_plot(fiti, pred = growth_mindset, modx = schoolInstrumental, x.label ="Mindset of Music Ability", y.label ="Academic Music Self-Concept", legend.main ="School Instrumental Ensemble")
Warning: 0.471724025973107 is outside the observed range of schoolInstrumental
Code
interact_plot(fiti, pred = schoolInstrumental, modx = growth_mindset, x.label ="School Instrumental", y.label ="Academic Music Self-Concept", legend.main ="Mindset of Music Ability")
Visual inspection of the data suggests coefficients for mindset of music ability, active engagement, school instrumental ensemble, extracurricular school music club, private lesson, and self-teaching may not be significantly different. Hence, I investigated further by conducting Wald Tests between each pair of predictors.
Linear hypothesis test:
privateLesson - selfTaught = 0
Model 1: restricted model
Model 2: musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental +
schoolExtracurricular + privateLesson + selfTaught + grade +
school + gender
Res.Df RSS Df Sum of Sq F Pr(>F)
1 490 2865.0
2 489 2815.1 1 49.905 8.6688 0.003391 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Wald tests results indicated that the difference in coefficients between mindset of music ability and extracurricular school music club, mindset of music ability and private lesson, active engagement and extracurricular school music club, active engagement and self-teaching, school instrumental ensemble and self-teaching, and private lesson and self-teaching were not statistically significant, suggesting those predictors had the same effect on students’ global music self-concept. Coefficients for all other pairs were statistically different from each other.
Joint Significance (Wald Test)
Code
wt_ability_grade <- lmtest::waldtest(lm_ability, term ="grade")wt_ability_grade
#scatter plots with fitted regression linesplot1_ability <-ggplot(data, aes(x = growth_mindset_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Mindset of Music Ability",y ="Ability Music Self-Concept") +theme_ipsum()plot2_ability <-ggplot(data, aes(x = active_engagement_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Active Engagement",y ="Global Music Self-Concept") +theme_ipsum()plot3_ability <-ggplot(data, aes(x = selfTaught_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Self-Teaching",y ="Global Music Self-Concept") +theme_ipsum()plot4_ability <-ggplot(data, aes(x = schoolExtracurricular_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Extracurricular School Music Club",y ="Global Music Self-Concept") +theme_ipsum()plot5_ability <-ggplot(data, aes(x = schoolInstrumental_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="School Instrumental Ensemble",y ="Global Music Self-Concept") +theme_ipsum()plot6_ability <-ggplot(data, aes(x = privateLesson_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Private Lesson",y ="Global Music Self-Concept") +theme_ipsum()#Combine plots into a single display (grid layout)grid.arrange(plot1_ability, plot2_ability, plot3_ability, plot4_ability, plot5_ability, plot6_ability, ncol =2)
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
Alternative Scatter Plots Option
Code
#ALTERNATIVE TEST#scatter plots with fitted regression linesggplot(data, aes(x = growth_mindset_z, y = musical_ability_msc_z, color = school_music_elective)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(title ="Ability Music Self-Concept ~ Mindset of Music Ability",x ="Mindset of Music Ability",y ="Ability Music Self-Concept",color ="Enrolled in music elective?") +theme_ipsum()
`geom_smooth()` using formula = 'y ~ x'
Code
ggplot(data, aes(x = active_engagement_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Active Engagement",y ="Global Music Self-Concept") +theme_ipsum()
`geom_smooth()` using formula = 'y ~ x'
Code
ggplot(data, aes(x = selfTaught_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Self-Teaching",y ="Global Music Self-Concept") +theme_ipsum()
`geom_smooth()` using formula = 'y ~ x'
Code
ggplot(data, aes(x = schoolExtracurricular_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Extracurricular School Music Club",y ="Global Music Self-Concept") +theme_ipsum()
`geom_smooth()` using formula = 'y ~ x'
Code
plot5_ability <-ggplot(data, aes(x = schoolInstrumental_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="School Instrumental Ensemble",y ="Global Music Self-Concept") +theme_ipsum()
NON-ACADEMIC MSC - Stepwise Multiple Linear Regression
Intercept
Code
#Interceptlm_msc_no_ability_step_0 <-lm(formula = msc_no_ability ~1, data = data)lm_msc_no_ability_step_0
Call:
lm(formula = msc_no_ability ~ growth_mindset + active_engagement,
data = data)
Residuals:
Min 1Q Median 3Q Max
-28.6795 -5.5348 0.4917 5.7274 28.0049
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.28260 2.25993 9.860 < 2e-16 ***
growth_mindset 1.11784 0.15279 7.316 1.02e-12 ***
active_engagement 1.59886 0.07735 20.670 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.983 on 500 degrees of freedom
Multiple R-squared: 0.5496, Adjusted R-squared: 0.5478
F-statistic: 305.1 on 2 and 500 DF, p-value: < 2.2e-16
Code
# Calculate model fit improvementf_test_result_lm_msc_no_ability_1_2 <-anova(lm_msc_no_ability_step_1, lm_msc_no_ability_step_2)f_test_result_lm_msc_no_ability_1_2
Analysis of Variance Table
Model 1: msc_no_ability ~ growth_mindset
Model 2: msc_no_ability ~ growth_mindset + active_engagement
Res.Df RSS Df Sum of Sq F Pr(>F)
1 501 74815
2 500 40343 1 34472 427.24 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculate model fit improvementf_test_result_lm_msc_no_ability_2_3 <-anova(lm_msc_no_ability_step_2, lm_msc_no_ability_step_3)f_test_result_lm_msc_no_ability_2_3
Call:
lm(formula = msc_no_ability ~ growth_mindset + active_engagement +
schoolChoir + schoolInstrumental + schoolExtracurricular +
outsideChoir + outsideInstrumental + privateLesson + selfTaught +
grade + school + familyMusic + race + hispanic_latino + gender,
data = data)
Residuals:
Min 1Q Median 3Q Max
-27.8726 -5.4640 0.3922 5.7536 28.7744
Coefficients:
Estimate Standardized Std. Error t value
(Intercept) 24.4065909 NA 3.4535028 7.067
growth_mindset 1.1744169 0.2399478 0.1727753 6.797
active_engagement 1.4671971 0.5921259 0.0898091 16.337
schoolChoir -0.0398351 -0.0053960 0.2490325 -0.160
schoolInstrumental -0.3717361 -0.0449088 0.2787228 -1.334
schoolExtracurricular -0.1361955 -0.0152693 0.3121567 -0.436
outsideChoir 0.3230472 0.0230407 0.4661907 0.693
outsideInstrumental -0.1079797 -0.0068968 0.5314533 -0.203
privateLesson -0.3145132 -0.0470852 0.2293146 -1.372
selfTaught 0.7664320 0.1119425 0.2497520 3.069
grade7 -0.1127632 -0.0039222 1.0451684 -0.108
grade8 0.1171442 0.0042469 1.0440654 0.112
schoolB -0.2212530 -0.0082851 1.0264244 -0.216
schoolC -0.8822040 -0.0292941 1.1601239 -0.760
familyMusicYes -0.3262705 -0.0121840 0.8437192 -0.387
raceBlack or African American 1.9593949 0.0406179 2.6207337 0.748
raceSome Other Race 2.0222003 0.0296097 2.9907924 0.676
raceTwo or More Races 1.8872029 0.0364543 2.7166139 0.695
raceWhite 1.1917803 0.0376185 2.2274437 0.535
hispanic_latinoYes 0.5822449 0.0093004 2.0399304 0.285
genderMale -1.8684377 -0.0698732 0.8824830 -2.117
genderNon-Binary -0.0460634 -0.0003066 4.6226558 -0.010
genderPrefer not to say -3.1897261 -0.0299035 3.2713764 -0.975
Pr(>|t|)
(Intercept) 5.60e-12 ***
growth_mindset 3.17e-11 ***
active_engagement < 2e-16 ***
schoolChoir 0.87298
schoolInstrumental 0.18293
schoolExtracurricular 0.66281
outsideChoir 0.48868
outsideInstrumental 0.83908
privateLesson 0.17085
selfTaught 0.00227 **
grade7 0.91413
grade8 0.91071
schoolB 0.82942
schoolC 0.44737
familyMusicYes 0.69915
raceBlack or African American 0.45504
raceSome Other Race 0.49928
raceTwo or More Races 0.48759
raceWhite 0.59287
hispanic_latinoYes 0.77544
genderMale 0.03475 *
genderNon-Binary 0.99205
genderPrefer not to say 0.33003
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 8.961 on 480 degrees of freedom
Multiple R-squared: 0.5697, Adjusted R-squared: 0.55
F-statistic: 28.89 on 22 and 480 DF, p-value: < 2.2e-16
Code
# Calculate model fit improvementf_test_result_lm_msc_no_ability_3_4 <-anova(lm_msc_no_ability_step_3, lm_msc_no_ability_step_4)f_test_result_lm_msc_no_ability_3_4
#RESULTS TABLE Non-Academic MSC# Extract coefficients and statisticsresults_lm_msc_no_ability <-bind_rows(broom::tidy(lm.beta(lm_msc_no_ability_step_4), conf.int =TRUE))# Round output to two decimal pointsresults_lm_msc_no_ability$estimate <-round(results_lm_msc_no_ability$estimate,2)results_lm_msc_no_ability$std_estimate <-round(results_lm_msc_no_ability$std_estimate,2)results_lm_msc_no_ability$std.error <-round(results_lm_msc_no_ability$std.error,2)results_lm_msc_no_ability$statistic <-round(results_lm_msc_no_ability$statistic,2)results_lm_msc_no_ability$p.value <-round(results_lm_msc_no_ability$p.value,3)results_lm_msc_no_ability$conf.low <-round(results_lm_msc_no_ability$conf.low,2)results_lm_msc_no_ability$conf.high <-round(results_lm_msc_no_ability$conf.high,2)# Combine conf.low and conf.high into a single column called "95% CI" and change column namesresults_lm_msc_no_ability <- results_lm_msc_no_ability %>%mutate("95% CI"=paste0("[", conf.low, ", ", conf.high, "]")) %>%select(-conf.low, -conf.high) %>%# Remove the original conf.low and conf.high columnsrename(b = estimate, β = std_estimate,SE = std.error,t = statistic,p = p.value)# Change order of columnsresults_lm_msc_no_ability <- results_lm_msc_no_ability %>%select(term, b, SE, β, t, "95% CI", p)# Change output for p values undeer .001results_lm_msc_no_ability <- results_lm_msc_no_ability %>%mutate(p =ifelse(p <0.001, "<.001", p))# Change predictors' namesresults_lm_msc_no_ability <- results_lm_msc_no_ability %>%mutate(term =case_match(term, "(Intercept)"~"(Intercept)","growth_mindset"~"Growth mindset", "active_engagement"~"Active Engagement","schoolChoir"~"School choir","schoolInstrumental"~"School instrumental","schoolExtracurricular"~"School music club","outsideChoir"~"Choir (outside school)","outsideInstrumental"~"instrumental (outside school)","privateLesson"~"Private lesson","selfTaught"~"Self-teaching","grade7"~"Grade 7","grade8"~"Grade 8","schoolB"~"School B","schoolC"~"School C","familyMusicYes"~"Family music engagement","raceBlack or African American"~"Black or African American","raceSome Other Race"~"Other race","raceTwo or More Races"~"Two or more races","raceWhite"~"White","hispanic_latinoYes"~"Ethnicity","genderMale"~"Male","genderNon-Binary"~"Non-Binary","genderPrefer not to say"~"Gender not disclosed") )print(results_lm_msc_no_ability, n=100)
# A tibble: 23 × 7
term b SE β t `95% CI` p
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 (Intercept) 24.4 3.45 NA 7.07 [NA, NA] <.001
2 Growth mindset 1.17 0.17 0.24 6.8 [-0.1, 0.58] <.001
3 Active Engagement 1.47 0.09 0.59 16.3 [0.42, 0.77] <.001
4 School choir -0.04 0.25 -0.01 -0.16 [-0.49, 0.48] 0.873
5 School instrumental -0.37 0.28 -0.04 -1.33 [-0.59, 0.5] 0.183
6 School music club -0.14 0.31 -0.02 -0.44 [-0.63, 0.6] 0.663
7 Choir (outside school) 0.32 0.47 0.02 0.69 [-0.89, 0.94] 0.489
8 instrumental (outside school) -0.11 0.53 -0.01 -0.2 [-1.05, 1.04] 0.839
9 Private lesson -0.31 0.23 -0.05 -1.37 [-0.5, 0.4] 0.171
10 Self-teaching 0.77 0.25 0.11 3.07 [-0.38, 0.6] 0.002
11 Grade 7 -0.11 1.05 0 -0.11 [-2.06, 2.05] 0.914
12 Grade 8 0.12 1.04 0 0.11 [-2.05, 2.06] 0.911
13 School B -0.22 1.03 -0.01 -0.22 [-2.03, 2.01] 0.829
14 School C -0.88 1.16 -0.03 -0.76 [-2.31, 2.25] 0.447
15 Family music engagement -0.33 0.84 -0.01 -0.39 [-1.67, 1.65] 0.699
16 Black or African American 1.96 2.62 0.04 0.75 [-5.11, 5.19] 0.455
17 Other race 2.02 2.99 0.03 0.68 [-5.85, 5.91] 0.499
18 Two or more races 1.89 2.72 0.04 0.69 [-5.3, 5.37] 0.488
19 White 1.19 2.23 0.04 0.54 [-4.34, 4.41] 0.593
20 Ethnicity 0.58 2.04 0.01 0.29 [-4, 4.02] 0.775
21 Male -1.87 0.88 -0.07 -2.12 [-1.8, 1.66] 0.035
22 Non-Binary -0.05 4.62 0 -0.01 [-9.08, 9.08] 0.992
23 Gender not disclosed -3.19 3.27 -0.03 -0.98 [-6.46, 6.4] 0.33
Visual inspection of the data suggests coefficients for mindset of music ability, active engagement, and self-teaching may not be significantly different. Hence, I investigated further by conducting Wald Tests between each pair of predictors.
Linear hypothesis test:
active_engagement - selfTaught = 0
Model 1: restricted model
Model 2: msc_no_ability ~ growth_mindset + active_engagement + selfTaught +
gender
Res.Df RSS Df Sum of Sq F Pr(>F)
1 497 39843
2 496 39157 1 685.21 8.6795 0.00337 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Wald tests results indicated statistical difference between the coefficients of mindset of music ability and self-teaching, suggesting they had the same effect on students’ global music self-concept. Coefficients between the dyads mindset of music ability/active engagement and active engagement/self-teaching were statistically different from each other.
Joint Significance (Wald Test)
Code
wt_gender <- lmtest::waldtest(lm, term ="gender")print(wt_gender)
Wald test
Model 1: msc_no_ability ~ growth_mindset + active_engagement + selfTaught +
gender
Model 2: msc_no_ability ~ growth_mindset + active_engagement + selfTaught
Res.Df Df F Pr(>F)
1 496
2 499 -3 2.0287 0.109
Scatter Plots
Code
#scatter plots with fitted regression linesplot1_no_ability <-ggplot(data, aes(x = growth_mindset_z, y = msc_no_ability_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Mindset of Music Ability",y ="Non-Academic Music Self-Concept") +theme_ipsum()plot2_no_ability <-ggplot(data, aes(x = active_engagement_z, y = msc_no_ability_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Active Engagement",y ="Global Music Self-Concept") +theme_ipsum()plot3_no_ability <-ggplot(data, aes(x = selfTaught_z, y = msc_no_ability_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Self-Teaching",y ="Global Music Self-Concept") +theme_ipsum()#Combine plots into a single display (grid layout)grid.arrange(plot1_no_ability, plot2_no_ability, plot3_no_ability, ncol =2)
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
Source Code
---title: "Regression Analysis"author: Diego Pintoformat: htmlself-contained: truetoc: truecode-fold: truecode-tools: truetheme: light: cosmoeditor_options: chunk_output_type: inline---# Set Up```{r}# Set working directory pathsetwd("~/Desktop/R_projects/msc")# Clear Global Environmentrm(list =ls()) # Load packageslibrary(tidyverse)library(ggplot2)library(gridExtra)library(hrbrthemes)library(psych)library(corrr)library(lm.beta)library (jtools)library(broom.mixed)library(car)# Read .csv file into a tibbledata <-read_csv("data/middle_school_music_mindsets_clean_partial.csv")``````{r}#reclassify character variables as factorsreclassFactor <-c("grade","school","race","hispanic_latino","gender","currentSchoolMusic","familyMusic")data[reclassFactor] <-lapply(data[reclassFactor], as.factor)``````{r}# Reorder the levels for currentSchoolMusic and set "No music" as the reference leveldata$currentSchoolMusic <-relevel(data$currentSchoolMusic, ref ="No music")```# Preliminary Analyses## Multicollinearity - VIF (Variance Inflation Factors)I explored the predictors' variance using Variance Inflation Factors (VIF) to identify potential collinearity. I followed the same steps for the three models (three different outcome variables).### Global MSC as outcome - VIF (Variance Inflation Factors)```{r}#VIF - All 17 predictorscar::vif(lm(formula = msc_full ~ growth_mindset + active_engagement + musical_training + currentSchoolMusic + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))```------------------------------------------------------------------------Based on these results, musical_training, currentSchoolMusic, and schoolInstrumental should not be included in the model together due to their potential multicollinearity. To solve this issue, I removed one by one starting with the one with highest VIF, currentSchoolMusic.```{r}#VIF - After removing currentSchoolMusiccar::vif(lm(formula = msc_full ~ growth_mindset + active_engagement + musical_training + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))```------------------------------------------------------------------------After removing currentSchoolMusic, VIF values for musical_training and schoolInstrumental are lower, but musical_training still hows high.```{r}#VIF - After removing currentSchoolMusic and musical_training car::vif(lm(formula = msc_full ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))```------------------------------------------------------------------------After removing currentSchoolMusic and musical_training, there seems to be no collinearity between the predictors. Hence,I will not include those variables in the model.### Ability MSC as coutcome - VIF (Variance Inflation Factors)```{r}#VIF - All 17 predictorscar::vif(lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + musical_training + currentSchoolMusic + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))```------------------------------------------------------------------------Similar to the previous model, musical_training, currentSchoolMusic, and schoolInstrumental should not be included in the model together due to their potential multicollinearity. To solve this issue, I removed one by one starting with the one with highest VIF, currentSchoolMusic.```{r}#VIF - After removing currentSchoolMusiccar::vif(lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + musical_training + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))```------------------------------------------------------------------------After removing currentSchoolMusic, VIF values for musical_training and schoolInstrumental are lower, but musical_training still hows high.```{r}#VIF - After removing currentSchoolMusic and musical_training car::vif(lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))```------------------------------------------------------------------------After removing currentSchoolMusic and musical_training, there seems to be no collinearity between the predictors. Hence,I will not include those variables in the model.### Non-Academic MSC as Outcome - VIF (Variance Inflation Factors)```{r}#VIF - All 17 predictorscar::vif(lm(formula = msc_no_ability ~ growth_mindset + active_engagement + musical_training + currentSchoolMusic + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))```------------------------------------------------------------------------Consistent with the first two models, musical_training, currentSchoolMusic, and schoolInstrumental should not be included in the model together due to their potential multicollinearity. To solve this issue, I removed one by one starting with the one with highest VIF, currentSchoolMusic.```{r}#VIF - After removing currentSchoolMusiccar::vif(lm(formula = msc_no_ability ~ growth_mindset + active_engagement + musical_training + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))```------------------------------------------------------------------------After removing currentSchoolMusic, VIF values for musical_training and schoolInstrumental are lower, but musical_training still hows high.```{r}#VIF - After removing currentSchoolMusic and musical_training car::vif(lm(formula = msc_no_ability ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data))```------------------------------------------------------------------------After removing currentSchoolMusic and musical_training, there seem to be no collinearity between the predictors. Hence,I will not include those variables in the model.{{< pagebreak >}}# AssumptionsBefore building the linear regression models, I examined assumptions of linearity, homoscedasticity, and other assumptions. I also examine if using the music self-concept composite score as an outcome would yield different results compared to the musical ability/music self-concept subscale.```{r}#Global MSC Final Modelfinal_model_msc_full <-lm(formula = msc_full ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data)#Ability MSC Final Modelfinal_model_msc_ability <-lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data)#No Ability MSC Final Modelfinal_model_msc_no_ability <-lm(formula = msc_no_ability ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data)```### Normality### Q-Q Plot```{r}qqPlot(final_model_msc_full)qqPlot(final_model_msc_ability)qqPlot(final_model_msc_no_ability)```### Kolmogorov-Smirnov```{r}# Global MSCks.test(residuals(lm(formula = msc_full ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data)), "pnorm")#Ability MSCks.test(residuals(lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data)), "pnorm")#No Ability MSCks.test(residuals(lm(formula = msc_no_ability ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data)), "pnorm")```### Shapiro-Wilk```{r}# Global MSCshapiro.test(residuals(lm(formula = msc_full ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data)))#Ability MSCshapiro.test(residuals(lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data)))#No Ability MSCshapiro.test(residuals(lm(formula = msc_no_ability ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data)))```Kolmogorov-Smirnov and Shapiro- Wilk tests showed that the data follows a nonnormal distribution. Individual examination of each variables' skewnewss and kurtosis identify three variables, but the distribution was not extreme.# Stepwise Multiple Linear RegressionI conducted a series of linear regression models to examine the relationships between multiple predictors and each of the three outcome variables, global MSC, ability MSC, and non-academic MSC. The aim was to address the research questions gradually, facilitating the interpretability of the findings by assessing the influence of successive categories of predictors on the models' improvement. I conducted the same steps and used the same predictors for each outcome variable.On Step 1, I introduced mindset of music ability as the sole predictor.To investigate the impact of various types of music engagement on music self-concept, I introduced two sets of variables in the following steps. Step 2 included active engagement, whereas Step 3 comprised variables related to duration of engagement in music activities in and out of school.In the following and final step (Step 4), I introduced background variables like grade, school, race, ethnicity, gender, and family music background.## GLOBAL MSC - Stepwise Multiple Linear Regression### Intercept```{r}#Interceptlm_msc_full_step_0 <-lm(formula = msc_full ~1, data = data)summary(lm_msc_full_step_0)```------------------------------------------------------------------------### Step 1```{r}#Step 1lm_msc_full_step_1 <-lm(formula = msc_full ~ growth_mindset, data = data)summary_lm_msc_full_step_1 <-summary(lm_msc_full_step_1)summary_lm_msc_full_step_1```------------------------------------------------------------------------Mindset of music ability significantly predicted global MSC, explaining 27% of its variance (adjusted *R^2^* = `r round(summary_lm_msc_full_step_1$adj.r.squared,4)`, *p* \< .01).### Step 2```{r}#Step 2lm_msc_full_step_2 <-lm(formula = msc_full ~ growth_mindset + active_engagement, data = data)summary_lm_msc_full_step_2 <-summary(lm_msc_full_step_2)summary_lm_msc_full_step_2# Calculate model fit improvementf_test_result_lm_msc_full_1_2 <-anova(lm_msc_full_step_1, lm_msc_full_step_2)f_test_result_lm_msc_full_1_2```------------------------------------------------------------------------In step 2, active engagement emerged as a significant predictor, collectively contributing to 62% of the variance in MSC (adjusted *R^2^* = `r round(summary_lm_msc_full_step_2$adj.r.squared,4)`, *p* \< .001 ). The addition of these variables significantly improved the model, *F* = `r round(f_test_result_lm_msc_full_1_2$F[2],2)`, *p* \< .001).### Step 3```{r}#Step 3lm_msc_full_step_3 <-lm(formula = msc_full ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught, data = data)summary_lm_msc_full_step_3 <-summary(lm_msc_full_step_3)summary_lm_msc_full_step_3# Calculate model fit improvementf_test_result_lm_msc_full_2_3 <-anova(lm_msc_full_step_2, lm_msc_full_step_3)f_test_result_lm_msc_full_2_3```------------------------------------------------------------------------The refined model explained 65% of the variance in music self-concept (adjusted *R^2^* = `r round(summary_lm_msc_full_step_3$adj.r.squared,4)`, p \< .01), representing a significant fit improvement, *F* = `r round(f_test_result_lm_msc_full_2_3$F[2],2)`, *p* \< .001).### Step 4```{r}#Step 4lm_msc_full_step_4 <-lm(formula = msc_full ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data)summary_lm_msc_full_step_4 <-summary(lm.beta(lm_msc_full_step_4))summary_lm_msc_full_step_4# Calculate model fit improvementf_test_result_lm_msc_full_3_4 <-anova(lm_msc_full_step_3, lm_msc_full_step_4)f_test_result_lm_msc_full_3_4```------------------------------------------------------------------------In the new model, gender was a significant predictor of global music self-concept, with male students differing significantly from female students (reference group). This extended model explained 68% of the variation in music self-concept (adjusted *R^2^* = `r round(summary_lm_msc_full_step_4$adj.r.square,4)`, *p* = .049), showing a significantly superior fit compared to the previous model (*F =* `r round(f_test_result_lm_msc_full_3_4$F[2],2)`, *p* \< .01).### Results Table```{r}#RESULTS TABLE Global MSC# Extract coefficients and statisticsresults_lm_msc_full <-bind_rows(broom::tidy(lm.beta(lm_msc_full_step_4), conf.int =TRUE))# Round output to two decimal pointsresults_lm_msc_full$estimate <-round(results_lm_msc_full$estimate,2)results_lm_msc_full$std_estimate <-round(results_lm_msc_full$std_estimate,3)results_lm_msc_full$std.error <-round(results_lm_msc_full$std.error,2)results_lm_msc_full$statistic <-round(results_lm_msc_full$statistic,2)results_lm_msc_full$p.value <-round(results_lm_msc_full$p.value,3)results_lm_msc_full$conf.low <-round(results_lm_msc_full$conf.low,2)results_lm_msc_full$conf.high <-round(results_lm_msc_full$conf.high,2)# Combine conf.low and conf.high into a single column called "95% CI" and change column namesresults_lm_msc_full <- results_lm_msc_full %>%mutate("95% CI"=paste0("[", conf.low, ", ", conf.high, "]")) %>%select(-conf.low, -conf.high) %>%# Remove the original conf.low and conf.high columnsrename(b = estimate, β = std_estimate,SE = std.error,t = statistic,p = p.value)# Change order of columnsresults_lm_msc_full <- results_lm_msc_full %>%select(term, b, SE, β, t, "95% CI", p)# Change output for p values undeer .001results_lm_msc_full <- results_lm_msc_full %>%mutate(p =ifelse(p <0.001, "<.001", p))# Change predictors' namesresults_lm_msc_full <- results_lm_msc_full %>%mutate(term =case_match(term, "(Intercept)"~"(Intercept)","growth_mindset"~"Growth mindset", "active_engagement"~"Active Engagement","schoolChoir"~"School choir","schoolInstrumental"~"School instrumental","schoolExtracurricular"~"School music club","outsideChoir"~"Choir (outside school)","outsideInstrumental"~"instrumental (outside school)","privateLesson"~"Private lesson","selfTaught"~"Self-teaching","grade7"~"Grade 7","grade8"~"Grade 8","schoolB"~"School B","schoolC"~"School C","familyMusicYes"~"Family music engagement","raceBlack or African American"~"Black or African American","raceSome Other Race"~"Other race","raceTwo or More Races"~"Two or more races","raceWhite"~"White","hispanic_latinoYes"~"Ethnicity","genderMale"~"Male","genderNon-Binary"~"Non-Binary","genderPrefer not to say"~"Gender not disclosed") )print(results_lm_msc_full, n=100)```------------------------------------------------------------------------### Overall Significance```{r}#Calculate overall significance/fitf_test_result_lm_msc_full_0_4 <-anova(lm_msc_full_step_0, lm_msc_full_step_4)f_test_result_lm_msc_full_0_4```------------------------------------------------------------------------The overall significance of the final model was confirmed through the F-test (*F =* `r round(f_test_result_lm_msc_full_0_4$F[2],2)`, *p* \< .01).### Post Hoc Analyses#### Equality of coefficients (Wald Test)Although several predictors showed significant effect on the outcome, their coefficient may or may not differ from each other significantly. To examine if the differences between the coefficients (as observed in values of *b* and β ) are significant, I visually inspected those coefficients.```{r}#Visually compare coefficientslm_full <-lm(formula = msc_full ~ growth_mindset + active_engagement + selfTaught + gender, data = data)jtools::plot_summs(lm_full)```------------------------------------------------------------------------Visual inspection of the data suggests coefficients for mindset of music ability, active engagement, and self-teaching may not be significantly different. Hence, I investigated further by conducting Wald Tests between each pair of predictors.```{r}wt_mindset_active <-linearHypothesis(lm_full, "growth_mindset = active_engagement ")wt_mindset_active```------------------------------------------------------------------------```{r}wt_mindset_selfTaught <-linearHypothesis(lm_full, "growth_mindset = selfTaught ")wt_mindset_selfTaught```------------------------------------------------------------------------```{r}wt_active_selfTaught <-linearHypothesis(lm_full, "active_engagement = selfTaught ")wt_active_selfTaught```------------------------------------------------------------------------Wald tests results indicated that the difference in coefficients between mindset of music ability and active engagement, mindset of music abilty and self-teaching, and active engagement and self-teaching were not statistically significant, suggesting those predictors had the same effect on students' global music self-concept.#### Joint Significance (Wald Test)```{r}wt_gender <- lmtest::waldtest(lm_full, term ="gender")print(wt_gender)```------------------------------------------------------------------------I employed the Wald test once more to delve deeper into the joint significance of the gender variable. The results showed a significant difference between the model incorporating gender and the one excluding it, *F* = `r wt_gender$F[2]`, *p* = \< .01. This finding indicates that gender stands as a significant predictor of music self-concept, underscoring its importance in the predictive model.### Scatter Plots```{r}#scatter plots with fitted regression linesplot1_full <-ggplot(data, aes(x = growth_mindset_z, y = msc_full_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Mindset of Music Ability",y ="Overall") +theme_ipsum()plot2_full <-ggplot(data, aes(x = active_engagement_z, y = msc_full_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Active Engagement",y ="Global Music Self-Concept") +theme_ipsum()plot3_full <-ggplot(data, aes(x = selfTaught_z, y = msc_full_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Self-Teaching",y ="Global Music Self-Concept") +theme_ipsum()plot4_full <-ggplot(data, aes(x = schoolExtracurricular_z, y = msc_full_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Extracurricular School Music Clubs",y ="Global Music Self-Concept") +theme_ipsum()#Combine plots into a single display (grid layout)grid.arrange(plot1_full, plot2_full, plot3_full, plot4_full, ncol =2)```------------------------------------------------------------------------## ACADEMIC MSC (Music Ability Only) - Stepwise Multiple Linear Regression### Intercept```{r}#Interceptlm_msc_ability_step_0 <-lm(formula = musical_ability_msc ~1, data = data)summary(lm_msc_ability_step_0)```------------------------------------------------------------------------### Step 1```{r}#Step 1lm_msc_ability_step_1 <-lm(formula = musical_ability_msc ~ growth_mindset, data = data)summary_lm_msc_ability_step_1 <-summary(lm_msc_ability_step_1)summary_lm_msc_ability_step_1```------------------------------------------------------------------------### Step 2```{r}#Step 2lm_msc_ability_step_2 <-lm(formula = musical_ability_msc ~ growth_mindset + active_engagement, data = data)summary_lm_msc_ability_step_2 <-summary(lm_msc_ability_step_2)summary_lm_msc_ability_step_2# Calculate model fit improvementf_test_result_lm_msc_ability_1_2 <-anova(lm_msc_ability_step_1, lm_msc_ability_step_2)f_test_result_lm_msc_ability_1_2```------------------------------------------------------------------------### Step 3```{r}#Step 3lm_msc_ability_step_3 <-lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught, data = data)summary_lm_msc_ability_step_3 <-summary(lm_msc_ability_step_3)summary_lm_msc_ability_step_3# Calculate model fit improvementf_test_result_lm_msc_ability_2_3 <-anova(lm_msc_ability_step_2, lm_msc_ability_step_3)f_test_result_lm_msc_ability_2_3```------------------------------------------------------------------------### Step 4```{r}#Step 4lm_msc_ability_step_4 <-lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data)summary_lm_msc_ability_step_4 <-summary(lm.beta(lm_msc_ability_step_4))summary_lm_msc_ability_step_4# Calculate model fit improvementf_test_result_lm_msc_ability_3_4 <-anova(lm_msc_ability_step_3, lm_msc_ability_step_4)f_test_result_lm_msc_ability_3_4```------------------------------------------------------------------------### Results Table```{r}#RESULTS TABLE: Ability MSC# Extract coefficients and statisticsresults_lm_msc_ability <-bind_rows(broom::tidy(lm.beta(lm_msc_ability_step_4), conf.int =TRUE))# Round output to two decimal pointsresults_lm_msc_ability$estimate <-round(results_lm_msc_ability$estimate,2)results_lm_msc_ability$std_estimate <-round(results_lm_msc_ability$std_estimate,2)results_lm_msc_ability$std.error <-round(results_lm_msc_ability$std.error,2)results_lm_msc_ability$statistic <-round(results_lm_msc_ability$statistic,2)results_lm_msc_ability$p.value <-round(results_lm_msc_ability$p.value,3)results_lm_msc_ability$conf.low <-round(results_lm_msc_ability$conf.low,2)results_lm_msc_ability$conf.high <-round(results_lm_msc_ability$conf.high,2)# Combine conf.low and conf.high into a single column called "95% CI" and change column namesresults_lm_msc_ability <- results_lm_msc_ability %>%mutate("95% CI"=paste0("[", conf.low, ", ", conf.high, "]")) %>%select(-conf.low, -conf.high) %>%# Remove the original conf.low and conf.high columnsrename(b = estimate, β = std_estimate,SE = std.error,t = statistic,p = p.value)# Change order of columnsresults_lm_msc_ability <- results_lm_msc_ability %>%select(term, b, SE, β, t, "95% CI", p)# Change output for p values undeer .001results_lm_msc_ability <- results_lm_msc_ability %>%mutate(p =ifelse(p <0.001, "<.001", p))# Change predictors' namesresults_lm_msc_ability <- results_lm_msc_ability %>%mutate(term =case_match(term, "(Intercept)"~"(Intercept)","growth_mindset"~"Growth mindset", "active_engagement"~"Active Engagement","schoolChoir"~"School choir","schoolInstrumental"~"School instrumental","schoolExtracurricular"~"School music club","outsideChoir"~"Choir (outside school)","outsideInstrumental"~"instrumental (outside school)","privateLesson"~"Private lesson","selfTaught"~"Self-teaching","grade7"~"Grade 7","grade8"~"Grade 8","schoolB"~"School B","schoolC"~"School C","familyMusicYes"~"Family music engagement","raceBlack or African American"~"Black or African American","raceSome Other Race"~"Other race","raceTwo or More Races"~"Two or more races","raceWhite"~"White","hispanic_latinoYes"~"Ethnicity","genderMale"~"Male","genderNon-Binary"~"Non-Binary","genderPrefer not to say"~"Gender not disclosed") )print(results_lm_msc_ability, n=100)```------------------------------------------------------------------------### Overall Significance```{r}#overall significance/fitf_test_result_lm_msc_ability_0_4 <-anova(lm_msc_ability_step_0, lm_msc_ability_step_4)f_test_result_lm_msc_ability_0_4```------------------------------------------------------------------------### Interaction (Growth Mindset X Skills-Related Music Engagement Activities)```{r}#Add interactions calculation to lm_msc_ability_step_4lm_msc_ability_interactions <-lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender + growth_mindset:schoolChoir +#interactions growth_mindset:schoolInstrumental + growth_mindset:schoolExtracurricular + growth_mindset:outsideChoir + growth_mindset:outsideInstrumental + growth_mindset:privateLesson + growth_mindset:selfTaught, data = data)summary_lm_msc_ability_interactions <-summary(lm.beta(lm_msc_ability_interactions))summary(lm_msc_ability_interactions)# Calculate model fit improvementf_test_result_lm_msc_ability_4_interactions <-anova(lm_msc_ability_step_4, lm_msc_ability_interactions)f_test_result_lm_msc_ability_4_interactions``````{r}#Edit previous model to only include significant interactions (growth_mindset*schoolInstrumental)lm_msc_ability_interactions_short <-lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender + growth_mindset:schoolInstrumental, data = data)summary_lm_msc_ability_interactions_short <-summary(lm.beta(lm_msc_ability_interactions_short))summary(lm_msc_ability_interactions_short)# Calculate model fit improvementf_test_result_lm_msc_ability_4_interactions_short <-anova(lm_msc_ability_step_4, lm_msc_ability_interactions_short)f_test_result_lm_msc_ability_4_interactions_short``````{r}#Interaction plot, Source: https://cran.r-project.org/web/packages/interactions/vignettes/interactions.htmllibrary(interactions)library(jtools) # for summ()#states <- as.data.frame(state.x77)fiti <-lm(musical_ability_msc ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender + growth_mindset:schoolInstrumental, data = data)summ(fiti)#Plotinteract_plot(fiti, pred = growth_mindset, modx = schoolInstrumental, x.label ="Mindset of Music Ability", y.label ="Academic Music Self-Concept", legend.main ="School Instrumental Ensemble")interact_plot(fiti, pred = schoolInstrumental, modx = growth_mindset, x.label ="School Instrumental", y.label ="Academic Music Self-Concept", legend.main ="Mindset of Music Ability")```### Post Hoc Analyses#### Equality of coefficients (Wald Test)```{r}#Visually compare coefficientslibrary (jtools)library(broom.mixed)lm_ability <-lm(formula = musical_ability_msc ~ growth_mindset + active_engagement + schoolInstrumental + schoolExtracurricular + privateLesson + selfTaught + grade + school + gender, data = data)jtools::plot_summs(lm_ability)```------------------------------------------------------------------------Visual inspection of the data suggests coefficients for mindset of music ability, active engagement, school instrumental ensemble, extracurricular school music club, private lesson, and self-teaching may not be significantly different. Hence, I investigated further by conducting Wald Tests between each pair of predictors.```{r}library(car)wt_ability_mindset_active <-linearHypothesis(lm_ability, "growth_mindset = active_engagement ")wt_ability_mindset_active```------------------------------------------------------------------------```{r}wt_ability_mindset_schoolInstrumental <-linearHypothesis(lm_ability, "growth_mindset = schoolInstrumental ")wt_ability_mindset_schoolInstrumental```------------------------------------------------------------------------```{r}wt_ability_mindset_schoolExtracurricular <-linearHypothesis(lm_ability, "growth_mindset = schoolExtracurricular ")wt_ability_mindset_schoolExtracurricular```------------------------------------------------------------------------```{r}wt_ability_mindset_privateLesson <-linearHypothesis(lm_ability, "growth_mindset = privateLesson ")wt_ability_mindset_privateLesson```------------------------------------------------------------------------```{r}wt_ability_mindset_selfTaught <-linearHypothesis(lm_ability, "growth_mindset = selfTaught ")wt_ability_mindset_selfTaught```------------------------------------------------------------------------```{r}wt_ability_active_schoolInstrumental <-linearHypothesis(lm_ability, "active_engagement = schoolInstrumental ")wt_ability_active_schoolInstrumental```------------------------------------------------------------------------```{r}wt_ability_active_schoolExtracurricular <-linearHypothesis(lm_ability, "active_engagement = schoolExtracurricular ")wt_ability_active_schoolExtracurricular```------------------------------------------------------------------------```{r}wt_ability_active_privateLesson <-linearHypothesis(lm_ability, "active_engagement = privateLesson ")wt_ability_active_privateLesson```------------------------------------------------------------------------```{r}wt_ability_active_selfTaught <-linearHypothesis(lm_ability, "active_engagement = selfTaught ")wt_ability_active_selfTaught```------------------------------------------------------------------------```{r}wt_ability_schoolInstrumental_schoolExtracurricular <-linearHypothesis(lm_ability, "schoolInstrumental = schoolExtracurricular ")wt_ability_schoolInstrumental_schoolExtracurricular```------------------------------------------------------------------------```{r}wt_ability_schoolInstrumental_privateLesson <-linearHypothesis(lm_ability, "schoolInstrumental = privateLesson ")wt_ability_schoolInstrumental_privateLesson```------------------------------------------------------------------------```{r}wt_ability_schoolInstrumental_selfTaught <-linearHypothesis(lm_ability, "schoolInstrumental = selfTaught ")wt_ability_schoolInstrumental_selfTaught```------------------------------------------------------------------------```{r}wt_ability_schoolExtra_privateLesson <-linearHypothesis(lm_ability, "schoolExtracurricular = privateLesson ")wt_ability_schoolExtra_privateLesson```------------------------------------------------------------------------```{r}wt_ability_schoolExtral_selfTaught <-linearHypothesis(lm_ability, "schoolExtracurricular = selfTaught ")wt_ability_schoolExtral_selfTaught```------------------------------------------------------------------------```{r}wt_ability_private_selfTaught <-linearHypothesis(lm_ability, "privateLesson = selfTaught ")wt_ability_private_selfTaught```------------------------------------------------------------------------Wald tests results indicated that the difference in coefficients between mindset of music ability and extracurricular school music club, mindset of music ability and private lesson, active engagement and extracurricular school music club, active engagement and self-teaching, school instrumental ensemble and self-teaching, and private lesson and self-teaching were not statistically significant, suggesting those predictors had the same effect on students' global music self-concept. Coefficients for all other pairs were statistically different from each other.#### Joint Significance (Wald Test)```{r}wt_ability_grade <- lmtest::waldtest(lm_ability, term ="grade")wt_ability_grade```------------------------------------------------------------------------```{r}wt_ability_school <- lmtest::waldtest(lm_ability, term ="school")wt_ability_school```------------------------------------------------------------------------```{r}wt_ability_gender <- lmtest::waldtest(lm_ability, term ="gender")wt_ability_gender```------------------------------------------------------------------------### Scatter Plots```{r }#scatter plots with fitted regression linesplot1_ability <- ggplot(data, aes(x = growth_mindset_z, y = musical_ability_msc_z)) + geom_point() + geom_smooth(method = "lm", se = TRUE) + labs(x = "Mindset of Music Ability", y = "Ability Music Self-Concept") + theme_ipsum()plot2_ability <- ggplot(data, aes(x = active_engagement_z, y = musical_ability_msc_z)) + geom_point() + geom_smooth(method = "lm", se = TRUE) + labs(x = "Active Engagement", y = "Global Music Self-Concept") + theme_ipsum()plot3_ability <- ggplot(data, aes(x = selfTaught_z, y = musical_ability_msc_z)) + geom_point() + geom_smooth(method = "lm", se = TRUE) + labs(x = "Self-Teaching", y = "Global Music Self-Concept") + theme_ipsum()plot4_ability <- ggplot(data, aes(x = schoolExtracurricular_z, y = musical_ability_msc_z)) + geom_point() + geom_smooth(method = "lm", se = TRUE) + labs(x = "Extracurricular School Music Club", y = "Global Music Self-Concept") + theme_ipsum()plot5_ability <- ggplot(data, aes(x = schoolInstrumental_z, y = musical_ability_msc_z)) + geom_point() + geom_smooth(method = "lm", se = TRUE) + labs(x = "School Instrumental Ensemble", y = "Global Music Self-Concept") + theme_ipsum()plot6_ability <- ggplot(data, aes(x = privateLesson_z, y = musical_ability_msc_z)) + geom_point() + geom_smooth(method = "lm", se = TRUE) + labs(x = "Private Lesson", y = "Global Music Self-Concept") + theme_ipsum()#Combine plots into a single display (grid layout)grid.arrange(plot1_ability, plot2_ability, plot3_ability, plot4_ability, plot5_ability, plot6_ability, ncol = 2)```### Alternative Scatter Plots Option```{r}#ALTERNATIVE TEST#scatter plots with fitted regression linesggplot(data, aes(x = growth_mindset_z, y = musical_ability_msc_z, color = school_music_elective)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(title ="Ability Music Self-Concept ~ Mindset of Music Ability",x ="Mindset of Music Ability",y ="Ability Music Self-Concept",color ="Enrolled in music elective?") +theme_ipsum()ggplot(data, aes(x = active_engagement_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Active Engagement",y ="Global Music Self-Concept") +theme_ipsum()ggplot(data, aes(x = selfTaught_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Self-Teaching",y ="Global Music Self-Concept") +theme_ipsum()ggplot(data, aes(x = schoolExtracurricular_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Extracurricular School Music Club",y ="Global Music Self-Concept") +theme_ipsum()plot5_ability <-ggplot(data, aes(x = schoolInstrumental_z, y = musical_ability_msc_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="School Instrumental Ensemble",y ="Global Music Self-Concept") +theme_ipsum()```------------------------------------------------------------------------## NON-ACADEMIC MSC - Stepwise Multiple Linear Regression### Intercept```{r}#Interceptlm_msc_no_ability_step_0 <-lm(formula = msc_no_ability ~1, data = data)lm_msc_no_ability_step_0```------------------------------------------------------------------------### Step 1```{r}#Step 1lm_msc_no_ability_step_1 <-lm(formula = msc_no_ability ~ growth_mindset, data = data)summary_lm_msc_no_ability_step_1 <-summary(lm_msc_no_ability_step_1)summary_lm_msc_no_ability_step_1```------------------------------------------------------------------------### Step 2```{r}#Step 2lm_msc_no_ability_step_2 <-lm(formula = msc_no_ability ~ growth_mindset + active_engagement, data = data)summary_lm_msc_no_ability_step_2 <-summary(lm_msc_no_ability_step_2)summary_lm_msc_no_ability_step_2# Calculate model fit improvementf_test_result_lm_msc_no_ability_1_2 <-anova(lm_msc_no_ability_step_1, lm_msc_no_ability_step_2)f_test_result_lm_msc_no_ability_1_2```------------------------------------------------------------------------### Step 3```{r}#Step 3lm_msc_no_ability_step_3 <-lm(formula = msc_no_ability ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught, data = data)summary_lm_msc_no_ability_step_3 <-summary(lm_msc_no_ability_step_3)summary_lm_msc_no_ability_step_3# Calculate model fit improvementf_test_result_lm_msc_no_ability_2_3 <-anova(lm_msc_no_ability_step_2, lm_msc_no_ability_step_3)f_test_result_lm_msc_no_ability_2_3```------------------------------------------------------------------------### Step 4```{r}#Step 4lm_msc_no_ability_step_4 <-lm(formula = msc_no_ability ~ growth_mindset + active_engagement + schoolChoir + schoolInstrumental + schoolExtracurricular + outsideChoir + outsideInstrumental + privateLesson + selfTaught + grade + school + familyMusic + race + hispanic_latino + gender, data = data)summary_lm_msc_no_ability_step_4 <-summary(lm.beta(lm_msc_no_ability_step_4))summary_lm_msc_no_ability_step_4# Calculate model fit improvementf_test_result_lm_msc_no_ability_3_4 <-anova(lm_msc_no_ability_step_3, lm_msc_no_ability_step_4)f_test_result_lm_msc_no_ability_3_4```------------------------------------------------------------------------### Results Table```{r}#RESULTS TABLE Non-Academic MSC# Extract coefficients and statisticsresults_lm_msc_no_ability <-bind_rows(broom::tidy(lm.beta(lm_msc_no_ability_step_4), conf.int =TRUE))# Round output to two decimal pointsresults_lm_msc_no_ability$estimate <-round(results_lm_msc_no_ability$estimate,2)results_lm_msc_no_ability$std_estimate <-round(results_lm_msc_no_ability$std_estimate,2)results_lm_msc_no_ability$std.error <-round(results_lm_msc_no_ability$std.error,2)results_lm_msc_no_ability$statistic <-round(results_lm_msc_no_ability$statistic,2)results_lm_msc_no_ability$p.value <-round(results_lm_msc_no_ability$p.value,3)results_lm_msc_no_ability$conf.low <-round(results_lm_msc_no_ability$conf.low,2)results_lm_msc_no_ability$conf.high <-round(results_lm_msc_no_ability$conf.high,2)# Combine conf.low and conf.high into a single column called "95% CI" and change column namesresults_lm_msc_no_ability <- results_lm_msc_no_ability %>%mutate("95% CI"=paste0("[", conf.low, ", ", conf.high, "]")) %>%select(-conf.low, -conf.high) %>%# Remove the original conf.low and conf.high columnsrename(b = estimate, β = std_estimate,SE = std.error,t = statistic,p = p.value)# Change order of columnsresults_lm_msc_no_ability <- results_lm_msc_no_ability %>%select(term, b, SE, β, t, "95% CI", p)# Change output for p values undeer .001results_lm_msc_no_ability <- results_lm_msc_no_ability %>%mutate(p =ifelse(p <0.001, "<.001", p))# Change predictors' namesresults_lm_msc_no_ability <- results_lm_msc_no_ability %>%mutate(term =case_match(term, "(Intercept)"~"(Intercept)","growth_mindset"~"Growth mindset", "active_engagement"~"Active Engagement","schoolChoir"~"School choir","schoolInstrumental"~"School instrumental","schoolExtracurricular"~"School music club","outsideChoir"~"Choir (outside school)","outsideInstrumental"~"instrumental (outside school)","privateLesson"~"Private lesson","selfTaught"~"Self-teaching","grade7"~"Grade 7","grade8"~"Grade 8","schoolB"~"School B","schoolC"~"School C","familyMusicYes"~"Family music engagement","raceBlack or African American"~"Black or African American","raceSome Other Race"~"Other race","raceTwo or More Races"~"Two or more races","raceWhite"~"White","hispanic_latinoYes"~"Ethnicity","genderMale"~"Male","genderNon-Binary"~"Non-Binary","genderPrefer not to say"~"Gender not disclosed") )print(results_lm_msc_no_ability, n=100)```------------------------------------------------------------------------### Overall Significance```{r}#overall significancef_test_result_lm_msc_no_ability_0_4 <-anova(lm_msc_no_ability_step_0, lm_msc_no_ability_step_4)print(f_test_result_lm_msc_no_ability_0_4) ```------------------------------------------------------------------------### Post Hoc Analyses#### Equality of coefficients (Wald Test)```{r}#Visually compare coefficientslibrary (jtools)library(broom.mixed)lm <-lm(formula = msc_no_ability ~ growth_mindset + active_engagement + selfTaught + gender, data = data)jtools::plot_summs(lm)```------------------------------------------------------------------------Visual inspection of the data suggests coefficients for mindset of music ability, active engagement, and self-teaching may not be significantly different. Hence, I investigated further by conducting Wald Tests between each pair of predictors.```{r}library(car)wt_mindset_active <-linearHypothesis(lm, "growth_mindset = active_engagement ")wt_mindset_active```------------------------------------------------------------------------```{r}wt_mindset_selfTaught <-linearHypothesis(lm, "growth_mindset = selfTaught ")wt_mindset_selfTaught```------------------------------------------------------------------------```{r}wt_active_selfTaught <-linearHypothesis(lm, "active_engagement = selfTaught ")wt_active_selfTaught```------------------------------------------------------------------------Wald tests results indicated statistical difference between the coefficients of mindset of music ability and self-teaching, suggesting they had the same effect on students' global music self-concept. Coefficients between the dyads mindset of music ability/active engagement and active engagement/self-teaching were statistically different from each other.#### Joint Significance (Wald Test)```{r}wt_gender <- lmtest::waldtest(lm, term ="gender")print(wt_gender)```------------------------------------------------------------------------### Scatter Plots```{r}#scatter plots with fitted regression linesplot1_no_ability <-ggplot(data, aes(x = growth_mindset_z, y = msc_no_ability_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Mindset of Music Ability",y ="Non-Academic Music Self-Concept") +theme_ipsum()plot2_no_ability <-ggplot(data, aes(x = active_engagement_z, y = msc_no_ability_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Active Engagement",y ="Global Music Self-Concept") +theme_ipsum()plot3_no_ability <-ggplot(data, aes(x = selfTaught_z, y = msc_no_ability_z)) +geom_point() +geom_smooth(method ="lm", se =TRUE) +labs(x ="Self-Teaching",y ="Global Music Self-Concept") +theme_ipsum()#Combine plots into a single display (grid layout)grid.arrange(plot1_no_ability, plot2_no_ability, plot3_no_ability, ncol =2)```------------------------------------------------------------------------
